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SUMMARY 


The status of an effort to increase the efficiency of calculating 
transient temperature fields in complex aerospace vehicle structures is 
described. The advantages and disadvantages of explicit and implicit 
algorithms are discussed. Explicit solution techniques require minimal 
computation per time step but have stability-limited step sizes. Implicit 
techniques permit larger step sizes because of better stability but 
require more computation per time step. A promising set of implicit 
algorithms, known as the GEAR package is described. Four test problems, 
used for evaluating and comparing various algorithms, have been selected 
and finite element models of the configurations are described. These 
problems include a Space Shuttle frame component, an insulated cylinder, a 
metallic panel for a thermal protection system and a model of the Space 
Shuttle Orbiter wing. Calculations were carried out using the SPAR finite 
element program, the MITAS lumped parameter program and a special purpose 
finite element program incorporating the GEAR algorithms. 

Results generally indicate a preference for implicit over explicit 
algorithms for solution of transient structural heat transfer problems 
when the governing equations are "stiff". Stiff equations are typical of 
many practical problems such as insulated metal structures and are 
characterized by widely differing time constants in the thermal response. 
In cases where implicit algorithms are appropriate, the GEAR algorithms 
offer high potential for providing increased computational efficiency. In 
some cases careful attention to modeling detail such as avoiding thin or 
short high-conducting elements can reduce the stiffness to the extent that 
explicit methods become advantageous. 


INTRODUCTION 

An effort is in progress at the NASA Langley Research Center to 
improve capability to predict and optimize the thermal -structural behavior 
or aerospace vehicle structures. The focus of this activity is on space 
transportation vehicles presently typified by the Space Shuttle Orbiter. 

A principal task is to significantly reduce the computing effort for 
obtaining transient temperature fields in the structure. This task is to 
be accomplished by incorporating the best state-of-the-art solution 
algorithms into general-purpose thermal analysis computer programs. 

Current activity is focused on evaluation and comparison of explicit and 
implicit solution algorithms. 

In reviewing current literature, a preference is evident amonj 
numerical analysis researchers for implicit algorithms for solution of 
stiff* sets of ordinary differential equations (ODE's). Many engineering 
analysts, however, prefer to use the longer-established explicit 


♦Stiff sets of ordinary differential equations are characterized by 
solutions with widely varying time-constants. The typical case is when 
the solution to the homogeneous problem has very small time constants 
compared to those of the forcing function (ref. 1). 



algorithms. A partial explanation for this dichotomy is that the full 
power cf the implicit approach has not been transferred from researchers 
to engineering analysts. 

In the explicit algorithms the set of temperatures at a given time is 
expressed as an explicit function of the set of previous temperatures in 
the structure. The time step (the difference between the present and 
previous times) is limited (often severely) in order that the technique be 
stable. In the implicit algorithms the present temperatures in the 
structure are interrelated through a set of algebraic equations (usually 
nonlinear) which are often costly to solve. For the commonly-used 
implicit algorithms there is no stability-imposed limitation on step 
size. The step size is limited by solution accuracy only, so that 
implicit algorithms can, in general, use much larger time steps than 
explicit algorithms. Because a single explicit time step is 
computationally faster than a single implicit time step the key to the 
advantageous use of implicit algorithms is to use the largest possible 
time step size. 

As presently implemented in thermal analysis computer programs, 
implicit algorithms generally require a user-specified fixed time step 
(refs. 2 to 6). The step size must be determined by trial, insight or 
other means. Because the user is usually unable to choose the largest 
possible time step at each time point the implicit algorithm is not used 
to maximum advantage. Furthermore, the solution must be repeated with a 
smaller time step in order to assess the error in the solution. The lack 
of automatic selection of step size based on a prescribed error tolerance 
has certainly delayed the full development cf the potential of implicit 
solution algorithms. 

The strategy being advocated in the solution of large problems by 
implicit methods is to have several alternate 'mplicit algorithms of 
varying order available and to automatically select both the largest 
possible time step ana the appropriate algorithm throughout the solution 
process (refs. 6,7). A promising set of algorithms, developed for the 
purpose of implementing the aforementioned strategy, is denoted the GEAR 
algorithms (refs. 7 to 10). Good 'performance of the GEAR algorithms has 
been demonstrated in applications to problems in structural dynamics, 
atmospheric pollution and hydrodynamics (ref. 7). These successes suggest 
the application of the GEAR techniques to transient thermal analysis. 

The purpose of the present paper is to describe the current status of 
ongoing evaluations and demonstrations of the use of explicit and implicit 
algorithms for transient thermal analysis of heated structures using the 
finite element method. A Shuttle frame test article, an insulated 
cylinder, a metallic multiwall thermal protection system panel, and a 
model of the Shuttle Orbiter wing are analyzed using the SPAR thermal 
analysis computer code (ref. 2). Comparisons between implicit and 
explicit algorithms are presented. The performance of the GEAR algorithms 
is evaluated for the cylinder problem. For benchmark checks the cylinder 
is also analyzed with the MITAS lumped parameter program (ref. 11). It is 
a characteristic of thermal analysis by finite element and lumped 
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parameter techniques that careful modeling can minimize stiffness in a 
problem and conversely, improper modeling can increase the stiffness. 
Since stiffness is one of the key factors in the performance of implicit 
and explicit algorithms, evaluations of these algorithms cannot be 
entirely separated from modeling considerations. Consequently the paper 
includes a limited study of the effects of modeling on the performance of 
the explicit and implicit algorithms. 

LIST OF SYMBOLS 


C capacitance matrix 

OT time step size 

e R error in numerical solution of the temperature at time t n 
truncation error of numerical integration method 
F right hand side of equation for transient problem* see Eq. (1) 

G right hand side of simplified transient problem, T = G(T,t) * C"V 

h n time step 

J Jacobian of system of differential equations = 9F/3T 
K conductivity matrix 

L length of a rod element 
Q thermal load vector 

q order of a multi step method 

R residual of the system of equations generated by the implicit method 

t time 

t n n-th time point 

T vector of temperatures 

Tj temperature at node i 

T 0 initial temperature at node 

a thermal diffusivity 

a coefficient in Adams-Moulton method, Eq. (20) 

B coefficient in backward difference method, Eq. (19) 

Superscripts 

i iteration number 

A dot represents differentiation with respect to time 


NATURE OF ALGORITHMS USED IN TRANSIENT THERMAL ANALYSIS 

A transient heat transfer problem when discretized by a finite 
element, finite difference or similar technique, is governed by the 
following system of equations 

CT = Q(T,t) -K(T,t)T = F(T,t) T(0) given (1) 

where F is generally a nonlinear function. It is usually impractical to 
obtain an analytical solution to Eq. (1) so that numerical integration 
methods are used. These methods obtain an approximation to the solution 
at discrete time points tj, t£» t 3 , . . . and are denoted time marching 
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schemes because the solution at a given time is obtained in terms of t.ie 
values at previous times. The simplest numerical integration technique is 
the Euler method which uses the first two terms in a Taylor series to 
predict T a', time t n+ j as 

T<W * T(t„) * h„ T (t„) 

- T(t n ) t h„ c- 1 F(T(t„),t n ) 

where 

= *n+l -t n (3) 

Euler's method is an example of an explicit integration technique, 
so-named because T(t n +i) is given explicitly in terms of known 
quantities. Another approach to the numerical integration of equation (1) 
is the backward difference method which is an example of an implicit 
method. In this approach 

T(t n +i) = T(t n ) + h n T(t n+1 ) 

= T(t n ) + h n C- 1 F(T(t n+1 ), t n+1 ) 

Equation (4) is a system of implicit equations for T(t n+ j), which is 
generally nonlinear and thus difficult to solve. The explicit algorithm 
is therefore easier to implement* and in general would be the best choice 
except for its stability properties. The term stability refers to the 
error propagation from one time step to the next. A method is unstable 
when an error in the solution at a time point is magnified at subsequent 
time points. 

To illustrate the problem of stability associated with explicit 
solution methods, it is instructive to examine the following simple 
example. Figure 1 shows a 2-node finite element where the temperature of 
node 2 is given as xt. Based on a linear temperature variation and lumpeJ 
capacitances, the differential equation for the temperature at node 1 is 

T = % (T 2 - T l ) (5) 

1 L 

where a is the diffusivity of the material and L is the length of the 
element. The exact solution to Eq. (5) is 

i 2 

Tj - -XL^/2a + Xt + (T q + XL^/2a)e ^ a ^ 1 ' (6) 


*The advantage of the explicit algorithm depends in part on the form of 
the capacitance matrix C. Usually, the capacitances are lumped and C is 
diagonal. In cases where C is not diagonal each step of the explicit 
method is much more costly. 
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which Is conposed of terms (first two terms) that vary slowly with respect 
to other terms (last). 

To study the stability of the explicit Euler method for this problem, 
assume an error e n in Ti(tJ and calculate the error e n+ j in Tj (t n +i) due 
to that error. From eqs. (2) and (5) 

^l(t n +l) s ^i(t n ) + h n Tj(t n ) 

2ah 

* ^i(t n ) [l-2h n «/L 2 ] + (-^) T 2 (t„) (7) 

From Eq. (7) 

e n+1 = (l-2h n a/L 2 )e n (8) 

For stability (that is, no error magnification) it is required that 

Ie n +il 4 I e n | (9) 

or 

h n 4 L 2 /« (10) 

For short or thin elements having high diffusivity, eq. (10) imposes a 
severe limit on the time step which can be taken by the explicit 
algorithm. For example, a 5mn thick aluminum element (a=7xl0“ 5 m 2 /s) 
requires 

h 4 (5 x 10“ 3 ) 2 /7 x 10-5 = 0.36 sec 

which is a very small time step when used for temperature histories of 
several hours. 

By contrast, the implicit integration method does not have a 
stability-limited time step. If the backward difference method, (eq. (4)), 
is used for eq. (5) one obtains 

T l( t n+l) * T l( £ n) * h n b (t n *l ) 

■ T,(t„) ♦ h„ (2«/L 2 )(T 2 (t„ +1 ) - TjIVj)) 

From Eq. (11), 

T 1 (t„+l ) ■ ET X (t„) ♦ h„ (2a/L 2 )T 2 (t ntl )]/(l + 2h„a/L 2 ) (12) . 

so that if the error in Ti(t n ) is e n , the error in Tj(t n+1 ) due to e n is 
e n+1 = e n /(l + 2h n a/L 2 ) (13) 
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From Eq. (13) it is clear that for any value of h n , e n +i will be smaller 
than e n . 

Another source of error, denoted the truncation error e t , is due to 
using only the first two terms in the Taylor series for estimating 
T(t n +i). It is easily shown that this part of the error for both the 
Euler method and the backward difference method is 

®t * ±!/ 2 "n 2 T(t„) (14) 

where the minus sign applies for Euler's method and the plus for the 
backward difference method. Since the exact solution to the example 
problem is known, the truncation error may be calculated exactly. Eqs. 

(6) and (14) lead to 

e t • ±(2a 2 h n 2 /L 4 )(T 0 + *L 2 / r 0 )e“ 2at/L2 (15) 

For small values of t the exponential is close to unity so that the 
following condition must be satisfied to avoid large errors. 

2 a 2 h n 2 /L 4 « 1 (16) 

or 

h n « L 2 /a (17) 

For large values of t, the exponential becomes very small and h can be 
large without causing a large ef In terms of Figure 1, small steps are 
required early in the temperature history but not later in the history. 
These conditions inmediately suggest the usefulness of variable time steps 
which are automatically selected according to the local behavior of the 
temperature response. 

This example problem exhibits the essential features of most 
transient conduction heat transfer problems with respect to the 
integration techniques, namely: 

(1) The thermal response may be divided into regions of slowly and 
rapidly varying temperatures. Steep transients accompany initial 
conditions or sudden changes in the heat load. 

(2) The rapidity of variation of the transient portion of the temperature 
history is proportional to the quantity t_ 2 /a. During such a 
transient, time steps much smaller than L z /a must be taken no matter 
what type of integration technique is used. 

(3) During a period of slowly-varying temperatures, large time steps may 
be taken by implicit integration techniques but explicit techniques 
must still use time steps which are less than L z /o. 

In mathematical terms, the sample problem is an example of a "stiff" 
problem. A stiff problem is one whose solution includes a slowly varying 
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function of time plus a transient function which changes rapidly. When 
explicit methods are applied to stiff problems, very small integration 
time steps must be taken even though the solution changes very slowly. 

For this reason stiff problems are usually best solved by implicit 
methods. The effort involved in solving a system such as Eq. (4) is 
usually cost-effective if a small number of large time steps are used. 

The Euler method and the backward difference methods are presented as 
representatives of a large class of explicit and implicit techniques, 
respectively. Higher-order methods typically use more previous 
information to predict the temperature at the current time and have 
truncation errors which are proportional to higher powers of h n . Such 
techniques are called multistep methods and their order is one less than 
the power of h n in the truncation error expression. The stability 
properties of multistep methods are similar to those of the Euler and 
backward difference methods. Most explicit methods are unstable for time 
steps much larger than LVo. Accordingly, thermal analysis computer 
programs generally select the explicit time step automatically based on 
the stability requirement. For implicit methods, accuracy primarily 
determines the step size, although stability may be a factor for highly 
nonlinear problems. Even for linear problems some implicit algorithms 
produce bounded oscillations if the time step size is too large (ref. 

14). Often, low-order implicit algorithms are less susceptible to these 
oscillations. It is concluded that a good package for integrating stiff 
systems of ordinary differential equations would be one which uses 
implicit methods and automatically selects the order and the step size 
based on desired accuracy. One package denoted the GEAR algorithms has 
these features and is discussed next. 


THE GEAR ALGORITHMS 

Several software packages based on the work of Gear have been 
developed for general use (ref. 7). The package most appropriate for 
application to finite element thermal analysis is denoted GEARIB*. This 
package is intended to solve systems of ordinary differential equations of 
the form 

C(T,t) T = F(T,t) (18) 

The package employs two classes of implicit multistep 
methods, Adams-Moulton and backward differences. For nonstiff equations 
the Adams-Moulton method of order one through twelve is used. This method 
has the general form 

T(t n +i) s T(t n ) + h n ^3 B-j T(tn+i_j) (19) 

i=0 


*An earlier and closely related software package denoted GEARB was 
developed to solve equations of the form T = G(T,t). In the present 
application G=C"^F. At this writing GEARIB has not been implemented and 
calculations have been performed using GEARB. 
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where q Is the order. For stiff equations the backward difference 
algorithms of orders one through five are used. These algorithms have the 
general form 

q 

T(t n +i) " h n e 0 T(t n +j) + ^ a-j T(t n +j_'j) (20) 

The coefficients and 3j are given in reference 10. The user 
selects the class of methods (Adams -Moulton or backward differences), and 
as described in reference 7 GEARIB automatically selects the appropriate 
time step and the order based on user specified error tolerance. It may 
seem surprising that implicit methods are used for both stiff and nonstiff 
problems. However, for nonstiff problems the set of algebraic equations 
associated with each time step may be solved very effectively and implicit 
methods often have smaller truncation errors than explicit methods of the 
same order. 

Use of the GEAR algorithms is explained by the backward difference 
algorithm (of order one). Applied to Eq. (18), eq (20) gives 

R = C[T(t n+ i) - T(t n )] - h n F(T(t n+ i), t n+ i) = 0 (21) 

This system of nonlinear algebraic equations is solved by the 
modified Newton's method. That is 

T 1 * 1 (t„ +1 ) -Tl(t„ +1 ) - [f^R (22) 

where 

-5T * c - V 


J =dF/^T is the Jacobian of the system at a previous time point and may 
be calculated according to one of four options specified by the user: 

Option 0 : The Jacobian is assumed to be the unit matrix. In this 

case the iteration represented by eq. (22) is very 
efficient to implement. However, it can be shown that it 
converges only for very small values of h n . The upper 
limit on h n is of the same order as that required for 
stability of an explicit method. As a result option 0 is 
similar in co^t and required step size to an explicit 
method. 

Option 1 : The Jacobian is calculated in a subroutine by the user. 

This is the recommended option for stiff problems. 

Option 2 : The same as option 1 except that the Jacobian is 

calculated by finite differences. This option is 
* intended for users who do not wish to supply a 
subroutine to calculate the Jacobian. 
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Option 3 : The Jacobian Is assumed to be diagonal and Is calculated 

by finite differences in GEAR. This option Is more 
costly per iteration than option 0 and more efficient 
than option 1 because aR/£T is diagonal. Its 
convergence properties are also in between options 0 and 
1 . 

For options 1, 2 and 3 the Jacobian is recalculated whenever the iterative 
solution of eq. (21) requires more than three iterations. The 
Adams-Moulton method has better accuracy but poorer stability than the 
backward difference method. Therefore for the stiff est problems backward 
differences with option 1 is generally reconmended and for nonstiff 
problems the Adams-Moulton method with option 0 is recommended (ref. 7). 

In this paper all applications of the GEAR algorithms use the backward 
difference method (eq. 20). 


DESCRIPTION OF TEST PROBLEMS AND RESULTS 
Insulated Shuttle Test Frame 

The first test problem used for algorithm evaluation is a Shuttle 
Orbiter frame analyzed and tested under transient heating as described in 
reference 12. The configuration shown in figure 2 consists of an aluminum 
frame surrounded by insulation. The principal purpose of the study of the 
configuration as discussed in reference 12, was to evaluate the thermal 
performance of the insulation during a simulated Shuttle flight. A 
secondary purpose was to evaluate the adequacy of thermal analysis 
procedures by analytical and test comparisons. 

The lumped parameter model received from the author of reference 12 
consists of a two Jimensional section of a symmetric half of the structure 
and contains 118 nodes (see figure 2b). The unknown temperatures are 
located at the centroids of the lumps. The lumped parameter model was 
converted to a finite element model for analysis using the SPAR program 
(ref. 2). The corresponding SPAR finite element model contains 149 grid 
points located at the ends or corners of the elements. The model contains 
148 elements including one-dimensional elements which account for 
conduction in the aluminum structure and radiation across the air gap and 
two-dimensional elements which model conduction in the insulation and 
across the gap. The difference in numbers of elements and grid points 
is due to the different modeling approaches of the two methods. 

Minor modifications were made to the finite element model following 
the conversion. These consisted of eliminating or consolidating some 
extremely thin or short finite elements in the aluminum structure in order 
to reduce the stiffness of the equations and to increase the allowable 
time step for the explicit solution algorithm. The properties of the 
aluminum structure are functions of temperature and the properties of the 
insulation are functions of temperature and pressure. Material properties 
are updated every 50 seconds. The pressure-dependence is treated in SPAR 
as time dependence since the pres sure -vs -time variation is known from the 
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trajectory data for the simulated flight conditions. The applied heating 
is specified by tabulations of temperatures at the outer surface of the 
insulation. 

The temperature history for the frame was ar' sing the SPAR 
explicit (Euler), and Implicit techniques - (Crz ' lichc " and backward 
differences). Comparisons of solution times are iver in ble 1. The 
explicit procedure using a time step of 0.16 s required 513 s of CPU* 
time. This time step was controlled by conduction chrough most of the 
aluminum elements along the center and front of the Fra' e. 

Solution time using the Crank -Nicholson algorithm varied from 380 s 
to 38 s as the time step was varied between 1.0 and 50 sec. The solution 
times for backward differences were close to those of Crank -Nicholson. As 
indicated in Table 1(b), there is very little loss of accuracy in either 
the structure or insulation temperatures with increased time step size. 

The conclusion is that there is over an order of magnitude difference in 
solution time between explicit and implicit solution techniques for the 
frame problem as modelled. 

T ? accuracy of the solutions by the various techniques is further 
assessed. Figure 3 contains temperature histories at a point in the outer 
layer of the aluminum structure corresponding to node 309 (see fig. 2b). 
The solid line in figure 3 represents the applied temperatures at the 
outer surface of the insulation (node 29). The dotted line shows 
temperatures obtained by the SPAR analysis. The SPAR temperatures are 
plotted as a single curve since there is little difference between the 
results. The dashed-dot line shows analytical results from the lumped 
parameter analysis of reference 12 which are also in close agreement with 
the SPAR temperatures. The circular symbols represent test data from 
reference 12. The closeness of all the results indicates that the models 
are adequate to simulate the temperature history in the test article. 

Insulated Cylinder 

Model Description. - For the next test problem, a configuration was sought 
which was larger (in terms of number of unkr.own temperatures) than the 
Shuttle frame and exhibited some of the characteristics of an insulated 
airframe structure. Also, a simple structure was sought so that a finite 
element model could be easily generated in a stand-alone program in which 
the GEARB algorithms could be incorporated and tested. These 
considerations led to the insulated aluminum cylindrical shell depicted in 
figure 4. The cylinder is 18m (720 in.) in length and 4.5m (180 in.) in 
diameter. The aluminum is 0.25cm (0.1 in.) thick and the insulation is 
5.0 cm (2.0 in) thick. The outer surface of the insulation is heated over 
a region which consists of one-third the length and half the 
circumference. The finite element model consists of a symmetric half of 
the cylinder and is composed of simple solid elements (K81 elements in 
SPAR). There are 39 elements along the cylinder length, 4 in the 


*A11 times are given for the Langley Research Center CYBER 173 computer 
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circumferential direction and 3 through the depth (2 elements in the 
insu'ation and one in the structure). Additionally, the outer surface of 
the insulation has quadrilateral elements (K41) which receive the heat 
load and quadrilateral radiation elements (R41) which radiate to space. 

As a result the model contains 800 grid points (hence unknown 
temperatures) and 650 elements. It is recognized that some features of 
the model are non optimum. For example modeling the thin aluminum layer 
with K81 elements and using large high-aspect-ratio elements are not 
considered good modeling practices. The effects of changing the model to 
assess these and other shortcomings are discussed later in the paper. The 
time-dependence of the heat load on the cylinder is plotted in figure 5. 
For all calculations material properties of the metal and insulation are 
temperature-dependent and are given in table 2. Material properties are 
updated every 200 seconds of the temperature history. 

Application of GEARB. - The GEARB algorithms were applied to this example 
using a special purpose finite element program which generates a finite 
element model of a cylinder using K81, K41, and R41 elements. The program 
contains the GEARB package and generates the matrix J and the vector G 
(see footnote, p. 7). Only the backward difference option is used because 
it is recommended for stiff sets of equations. The first set of 
calculations concerns selecting the best Jacobian option. Temperature 
histories in the cylinder were calculated for 2000 s using each of the 
four options with a specified relative error limit of 0.001. Solution 
times, integration step sizes, and the number of Jacobian evaluations are 
given in table 3(a). As expected for this stiff problem, the only useful 
options are the user-supplied Jacobian (Option 1) ara he finite 
difference Jacobian (Option 2). The degree of js is indicated by 

the small step size (0.045 s) required by the exp .-like option (Option 
0). In contrast, options 1 and 2 permit time s^p. of up to 93.4 s and 
average 25 s. 

To further investigate the ’acobian options, the problem was made 
less stiff by increasing the mecai thickness to 2.54cm (1.0 in) and the 
calculations were repeated. The results in Table 3(b) show that options 0 
and 3 are acceptable (in fact option 3 is better than option 2 for this 
case) but are not as effective as option 1. The results indicated in 
table 3 suggest that as less stiff problems are considered, options 0 and 
3 will become more effective. This is an important trend because use of 
options 0 and 3 requires less core storage than options 1 and 2 since a 
diagonal Jacobian is used in the latter options. 

To assess the effect of accuracy requirements on computation time and 
results, the thin cylinder was reanalyzed using option 1 and a relative 
accuracy of 0.01. The CPU time was reduced from 450 s to 263 s, the 
average time step increased to 69 s and the number of Jacobian evaluations 
reduced to 9. The largest difference in calculated temperatures resulting 
from the relaxed error tolerance was only 14K (out of 560K) for a point in 
the insulation. 

Application of SPAR .- The temperature history of the cylinder for 2000 s 
was computed with SPAR using th^. explicit Euler algorithm as well as the 
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Crank-Nicholson and backward difference implicit algorithms. Comparisons 
of solution times for the methods are shown in table 4. There was 
essentially no difference between solution times for Crank-Nicholson and 
backward differences and results are presented only for the former 
method. The explicit algorithm estimated a stability limit of 0.12 s, 
however, use of this time step gave an unstable solution. A time step v,, 
0.06 s was used successfully and the solution time was 12300 seconds. The 
time step was controlled by conduction through the thin aluminum 
structural elements. The value of LVa based on a value of L of 0.254 cm 
(0.1 in.) was 0.103 s. This example illustrates a case of a stiff problem 
made more stiff by a model which causes the explicit algorithm to use a 
small time step to compute the negligible temperature gradient through the 
aluminum skin. e implicit algorithm was used with time steps of 5, 10 
and 25 s and required solution times of 782, 569 and 507 s respectively. 
For a time step of 50 s the Implicit method failed to converge. The 
relatively small decrease in solution time between time steps of 10 s and 
25 s is noted and is due to two reasons. First, a major portion of the 
time is used in SPAR to regenerate the finite element matrices when 
material properties are temperature-dependent (see next section). Second, 
a larger time step often increases the number of iterations required to 
solve the implicit system (equation 21). 

Comparison between GEARB and Implicit Algorithms in SPAR. - Experience 
with the GEARB algorithms and those presently in SPAR plus comparisons of 
solution tines such as those in Table 4 suggests the following advantages 
of the GEARB methods: 

(i) the use of accuracy-controlled time steps frees the user from the 
need to determine time steps for achieving desired accuracy; 

(ii) The use of variable time steps permits much larger average time 
steps to be used; 

(iii) GEARB employs an efficient predictor (the algorithm that supplies 
the first guess to the solution of eq. 21) and therefore can save 
time by employing larger time steps. 

To gain additional insight into the benefits of variable time steps 
and order in GEARB, the cylinder was reanalyzed with the heat lead of 
figure 5 replaced by a step function having the same peak value c.s figure 
5. This load results in a rapid and highly nonlinear response during the 
first part of the temperature history. Temperature histories were 
computed using GEARB in the special purpose program and Crank-Nicholson in 
SPAR (with a time step of 25 s). Material properties were updated every 
50 seconds. The time step used in GEARB varied between 3.5 s a.id 308 s 
with an average time step of 47 s. The order of the algorithms varied 
between 3 and 1. The smaller time steps and higher orders were used early 
in the time history. Solution time for GEARB was 368 s. compared to 
1014 s for Crank-Nicholson. Additionally there was an error of 10K (out 
of 750K) in the Crank-Nicholson result at 50 s. 
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Effect of Modeling on Algorithm Performance. - The thickness of the 
cylinder used In the calculations Is deliberately chosen to be quite small 
(consistent with the Shuttle frame for example). For this thickness there 
Is no significant temperature gradient through the aluminum and there Is 
no need to use elements (e.g. K81) which account for the gradient. 
Additionally, It 1$ possible to replace the three-dimensional K81 elements 
In the Insulation with an assemblage of one dimensional conductors through 
the Insulation thickness. Two models that reflect these Idea? were 
generated. The first model (model II) replaced the 3-dimensional aluminum 
elements with 2 dimensional (K41) elements. The second model (model III) 
used the two-dimensional aluminum elements and one-dimensional insulation 
elements. Both models have a finer (3 elements instead of 2) 
representation through the insulation so as to preserve the total number 
of grid points at 800. The solution times with these models are given in 
the third and fourth columns of Table 5. They indicate that the changes 
in the model which reduce the stiffness enable the explicit algorithms to 
execute faster than the implicit algorithms. As noted in the table, the 
implicit algorithms in SPAR for models II and III did not converge for a 
time step of 25 s and for mode' I, the solution time was greater for 25 s 
than for 17 s. These are additional indications that the predictor used 
in conjunction with the iterative solution of eq.(21) may be deficient. 

Another aspect of the effect of modeling is comparison of results 
from finite element and lumped parameter models. For this purpose, the 
MITAS lumped parameter computer program (ref. 11) was applied to the 
analysis of the cylinder. The finite element model I was converted to a 
lumped parameter model by use of the CINGEN program (ref. 13).* The 
resulting lumped parameter model contained 625 nodes as compared to 800 
grid points in the finite element model. Recall the unknown MITAS 
temperatures are located only at the centroids of each lump or element. 
Temperature histories were obtained using MITAS with the explicit (forward 
differences) and implicit (Crank -Nicholson and backward differences) 
methods. 

MITAS computation times are shown in the last column of table 5. 
Because none of the SPAR models is equivalent to the MITAS model in terms 
of the number of unknown temperature or nodal connections, no direct 
comparison of MITAS and SPAR solution times are appropriate. However, 
some trends evident in table 5 are noted. The MITAS model is not 
particularly stiff as evidenced by the large time step used in the 
explicit solution technique. The modified SPAR models which begin to 
resemble the MITAS model in certain respects are also less stiff and favor 
explicit algorithms. Of particular importance is the decrease in solution 
time of each program due to increased step size. Specifically note in 
trble 5 the large improvement in solution time between a time step of 10 
and 25 seconds in MITAS compared to the much smaller decrease in SPAR. 


♦CINGEN did not properly account for two-material conductors. These had 
to be input manually. 
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This is primarily due to the fact that SPAR regenerates the element 
conductivity matrices each time the temperature-dependent conductivities 
are updated. This results in high computation time especially for the 
solid (K81) elements. An alternative which can be easily implemented for 
isotropic elements (and is equivalent to what is done in MITAS) is to 
generate the matrices once and multiply each matrix by the latest updated 
conductivity. Presently the time used for the matrix regeneration in SPAR 
tends to mask some of the benefits of using implicit methods. Namely for 
larger time steps, the matrix regeneration time becomes a predominant 
portion of the total solution time. 

Figure 6 contains temperature histories of a point in the cylinder 
computed by the implicit and explicit techniques for all three SPAR 
models, plus GEARB and MITAS. Model II is considered to be best of the 
models being compared and thus the temperatures represented by the dotted 
line are thought to be the most accurate. These results are bracketed by 
results from model I and MITAS (from above) and by model II (from below). 
There are negligible differences between temperatures from the implicit 
and explicit solutions for any given model. Also GEARB produce" the same 
results as SPAR for model I. Results from model II and III are different 
from that of model I because of the extra layer of elements through the 
insulation. The MITAS tenperature history agrees well with that of model 
I (on which the MITAS model is based) except for some differences 
beginning at 1400 s. 

Multi wall Thermal Protection System Panel 

The next example problem is one which grew out of a study of the 
thermal performance of a titanium multiwall thermal protection system 
(TPS) panel which is under study for future use on space transportation 
systems (ref. 15). The configuration as depicted in figure 7(a), consists 
of alternating layers of flat and dimpled sheets fused at the crests to 
form a sandwich. The representation of a typical dimpled sheet is shown 
in figure 7(b). For the purpose of this analysis, it is assumed that the 
heat load does not vary in directions parallel to the plane of the panel. 
This assumption in addition to the regular geometry of the structure leads 
to the modeling simplification wherein only a triangular prismatic section 
of the panel needs to be modeled (fig 7(a)). The intersection of this 
prism with a typical dimpled layer is indicated by the shaded triangle in 
figure 7(b). 

The finite element model shown in figure 8 contains 333 grid points 
located on nine titanium sheets (5 horizontal and 4 inclined). The model 
contains 283 triangular and quadrilateral metal conduction elements, 264 
solid air conduction elements which account for gas conduction between the 
layers and 544 triangular and quadrilatral radiation elements which 
account for radiation heat transfer between adjacent horizontal and 
inclined sheets. Thermal properties of titanium and air are functions of 
temperature. Radiation exchange (view) factors were computed and supplied 
to SPAR using the TRASYS II computer program (ref. 16). 
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The temperature history of the panel In response to an imposed 
transient temperature at the outer surface of the panel was computed for 
2000 s. Results were obtained with SPAR using explicit. Crank 'Nicholson 
and backward difference algorithms. Solution-time comparisons are 
presented In table 6. The explicit algorithm required a time step of 
0.03 s. This time step was dictated by conduction of heat through the 
short heat paths between the verticies of adjacent triangular layers and 
indicates that this is an extremely stiff problem. Required solution time 
for the explicit algorithm was estimated to be 42000 s.* 

The Crank-Nicholson solution was carried out using time steps of 1 
and 5 s which led to solution times of 1811 and 675 s respectively. 
Backward differences was used with the same time steps and had solution 
times of 1772 and 703 s respectively. This example shows most 
dramatically the potential advantages of using implicit algorithms for 
thermal analysis of stiff problems. A plot of typical temperature 
histories for a point midway through the panel and the primary structure 
are shown in figure 9 along with the applied outer surface temperature. 

The results were obtained by the implicit algorithm with a time step of 5 
s and are identical to results using a time step of 1 s. 

Shuttle Wing 

The last example problem is a model of the Space Shuttle Orbiter 
wing. The model shown in figure 10 is based on a coarse (418 grid point) 
model and augmented by insulation attached to the upper and lower 
surfaces. The structure is modeled by rod, triangular and quadrilateral 
elements (K21,K31,K41 in SPAR terminology). The external insulation on 
each surface is modeled by five layers of solid triangular prismatic (K61 ) 
elements. The complete model contains 2508 grid points, 1400 one-and 
two-dimensional elements in the structure and 2700 solid elements in the 
insulation. Thermal properties of the alumini*m structure are 
temperature-dependenc; thermal properties of the insulation are 
temperature and time-dependent. 

For the purpose of this analysis, the applied heating on the wing is 
represented by a time-dependent temperature applied to the external 
surface of the insulation on the under side of the wing. The shape of 
this curve shown as the solid line in figure 11 is roughly indicative of 
atmospheric reentry heating. The temperature history of the wing for 4500 
seconds was computed using the SPAR explicit algorithm. 
Temperature-dependent properties were updated every 100 seconds of the 
temperature history. For this problem the explicit algorithm was able to 
use a large time step of 100 s. (The 100 s time step was dictated by the 
need to periodically update temperature-dependent material properties and 
not by stability requirements.) The time step is due to the coarse 
modeling of the structure which did not include the thin, high-conducting 
or radiating elements present in the previous models. Figure 11 shows the 


*To conserve computer resources the solution was terminated after 400 s of 
the temperature history for which 8400 CPU s were required. 
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temperature histories of a point on the structure and a point In the 
Insulation 1/5 of the distance through the Insulation at a typical cross 
section through the wing. The solution time for the explicit algorithm 
was 8600 s. Next the Implicit (Crank-NIcholson) solution algorithm was 
applied to the wing using a step size of 100 s. Over 5000 seconds of CPU 
time were used without completing the first time step. It was determined 
that the excessive slowness was due to the poor banding of the matrix 
equations which are solved as part of the Implicit technique (represented 
by eq. 21). The grid point decomposition sequence was changed In such a 
way as to greatly reduce the band width. The implicit solution was 
repp 'ted with the result that three time steps were completed using 1100 
seconds of CPU time. The solution was terminated after this point to 
conserve coiqputer resources. Extrapolating these values gives an estimate 
of 16500 CPU seconds to complete the 4500-second temperature history of 
the winy. Thus the implicit algorithm requires about twice the solution 
time as the explicit algorithm. This wing problem is a case where because 
of low stiffness the explicit algorithm is the best choice. It also shows 
that when using implicit methods, the analyst should be aware of the 
importance of proper banding of the matrices and careful grid point 
numbering. 


Choice of Explicit or Implicit Algorithms 

Two main factors determine whether explicit or implicit algorithms 
are more effective for solving a structural heat transfer problem. These 
are (1) stiffness of the system and (2) the connectivity of the model. 
These are now discussed in detail. 

Stiffness of the ODE system. -The stiffer the ODE system is, the more 
likely it is that the implicit algori* 1 s will be more efficient than 
explicit algorithms. In many cases ca ful and judicious modeling of the 
thermal problem can reduce the stiffness of the resulting system. 

However, the use of implicit algorithms can help the analyst avoid the 
added effort required tor such a judicious and careful modelling. 

The stiffness -yf the system depends primarily on the smallest time 
constant (LV<*) of the elements. Radiation and convection effects 
increase the stiffness of the model because they increase the conductance 
without affecting the capacitance, she stiffness of the system also 
depends on the applied heat loads. The system is stiff if these loads 
change much more slowly than the smallest time constant of the model. If 
the loads change very rapidly, small time steps are required to follow the 
response for both explicit and implicit algorithms. The explicit 
algorithms most likely will be more efficient in this case. 

Effect of connectivity of the model. -The disadvantage of an implicit 
method is associated with the need to solve a (generally nonlinear) system 
of equations such as eq. (21) at every time step. The use of the modified 
Newton method converts this problem to one of solving a series of linear 
systems. In SPAR, the linear system is solved by Gaussian elimination and 
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In MITAS by an iterative method. When the system of equations is poorly 
banded, Gaussian elimination is an ineffective solution method. 

Therefore, In SPAR it may be anticipated that the implicit methods become 
less attractive for problems which are poorly banded due to poor nodal 
numbering, or the inherent properties of the model. Problems with 
inter-element radiation, for example, tend to be poorly banded because 
non-adjacent grid points are coupled by radiation. 

The insulated cylinder problem (model I) is used to demonstrate the 
effect of banding on the implicit algorithms. The problem was originally 
modeled with 3 elements through the thickness, 4 in the circumferential 
direction and 39 in the axial direction. This model has 800 nodes with a 
band width of 51. The cylinder was remodeled with 9 elements through the 
thickness, 7 in the circumferential direction and 9 in the axial 
direction. This model also has 800 nodes but the band width is increased 
to 182. The solution time using the implicit algorithms for a time step 
of 5 seconds was 2670 seconds as compared to 782 seconds for the narrow 
band width cylinder (see table 5). This suggests that an inplicit scheme 
may certainly become less efficient if Gaussian elimination is used as the 
solution strategy and the system is poorly banded. 


CONCLUDING REMARKS 

This paper discusses the status of an effort to obtain increased 
efficiency in calculating transient temperature fields in complex 
aerospace vehicle structures. Explicit solution techniques which require 
minimal computation per time step and implicit techniques which permit 
larger time steps because of better stability are reviewed. A promising 
set of implicit solution algorithms, known as the GEARB and GEARIB 
packages are described. Four test problems for evaluating the algorithms 
have been selected and finite element models of each one are described. 

The problems include a Shuttle frame component, an insulated cylinder, a 
metallic panel for a thermal protection system and a model of the Space 
Shuttle Orbiter wing. Calculations were carried out using the SPAR finite 
element program, a special purpose finite element program incorporating 
the GEARB algorithms, and for checking purposes the MITAS lumped parameter 
program. 

Results generally indicate that implicit algorithms are more efficient 
than explicit algorithms for solution of transient structural heat 
transfer problems when the governing equations are stiff. Stiff equations 
are typical of many practical problems such as insulated metal structures 
and are characterized by widely differing time constants and cause 
explicit methods to take very small time steps. In those cases where 
implicit algorithms are appropriate, the GEARB and GEARIB algorithms offer 
high potential for obtaining the increased computational efficiency. 

Studies were also made of the effect on algorithm performance of 
different models of the same cylinder test problem. These studies 
revealed that the stiffness of the problem is highly sensitive to modeling 


17 



details and that careful modeling can reduce the stiffness of the 
resulting equations to the extent that explicit methods are advantageous. 
Since Implicit algorithms are less Influenced by stiff ness -related 
modeling details, use of these algorithms can save the analyst a certain 
amount of model refinement effort. Finally, wide-banding of the matrix 
equations of the finite element model either due to non-optlmal grid-point 
numbering or high connectivity (due for example to radiation) may decrease 
the advantage of implicit methods. 
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Table 1. - Performance of Various Algorithms for Transient Thermal 
Analysis of Shuttle Frame 

(a) Solution Time Comparison 


EXPLICIT 

IMPLICIT 

Euler 

Crank -Nicholson 

Backward 

Difference 

Time Step 

Solution 

Time Step 

Solution 

Time Step 

Solution 

(s) 

Time (s) 

(s) 

Time (s) 

(s) 

Time (s) 

0.16 

513 

1 

380 

1 

357 



10 

65 

10 

68 



25 

48 

25 

49 



50 

38 

50 

41 


(b) ffect of Time Step on Accuracy 
of Implicit Algorithms 


Step Si ze 
(s) 

Temp, of Node 309** at 1200 s 

Temp, of Node 409** at 1200 s 

K 

°F 

K 

°F 

1.0 

436.2 

325.2 

528.4 

mm 

10.0 

436.2 

325. 1 

528.4 


25.0 

436.1 

325.0 

528.3 

mm 

50.0 

437.2 

437.0 

528.6 


0.16* 

437.8 

328.0 

529.0 

IBiSfll 


* Explicit Algorithm 

** See figure 2(b) 
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Table 2. - Material Properties for Insulated Cylinder 

(a) Insulation: p = 160 kg/m^ (.00582 lbm/in^) 


T 

C 

k 

K 

°R 

J/kg-°C 

Btu/lbm-°R 

w/m-°C 

Btu/in-s-°R 

456 

360 

523 

0.125 

.0381 

5. 1x10-7 

622 

660 





.0546 

7.3 

733 

860 





.0711 

9.5 

844 

1060 





.0898 

1.2x10-6 

956 

1260 





.112 

1.5 

1067 

1460 





.142 

1.9 

1778 

1660 





.180 

2.4 


(b) Aluminum: p = 2770 kg/m 

3 (.101 Ibrn/in^) 


456 

360 

769 

0.184 

99.5 

.00133 

557 

560 

861 

.206 

125.0 

.00167 

622 

660 

903 

.216 

138.0 

.00185 

678 

760 

937 

.224 

154.8 

.00207 

733 

860 

974 

.233 

171.3 

.00229 

739 

960 

1012 

.242 

178.8 

.00239 

844 

1060 

1045 

.250 

181.1 

.00242 
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Table 3. - Effect on Soluction Time of Various GEARB 

Options for Jacobian Evaluation for Insulated 
Aluminum Cylinder 

(a) 0.254 cm (0.1 In.) Aluminum thickness 


Jacobian 

Option 

Unit 

Matrix 

(0) 

User 

Supplied 

(1) 

Finite 

Difference 

(2) 

Finite 
Diff (Diag) 
(3) 

CPU Time* 
(*) 

30,000** 

45C 

955 

10,000*** 

Step Size- 

Range 

Average 

0.045 

15-93.4 

25.0 

15.0-93.4 

25.0 

0.8 

Number of 

Jacobian 

Evaluations 

0 

17 

17 

2600 


i 1 ■ ■ ■' 1 1 1 ■ .i . , . . - , . _ , _ . i 

(b) 2.54 cm (1.0 in.) Aluminum thicknej 

>s 

CPU Time* 
(s) 

1075 

i 

402 

810 

703 

Step Size- 

Range 

Average 

1.6-14.8 

2.2 

13.8-83.2 

30.0 

13.8-83.2 

30.8 

2.2-21.9 

8.1 

Number of 

Jacobian 

Evaluations 

0 

14 

14 

223 


* For CDC CYBER 173 Computer 

** Estimate based on 1540 s for 100 s of History 

*** Estimate based on 1028 s for 200 s of History 
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Table 4. - Solution Times for Various Algorithms 
In Transient Thermal Analysis of 
Insulated Cylinder 


Explicit 

Euler 

implicit 

Crank-Nichol son/ 
Backward Difference 

GEARB 

Time Step 

Solution 

Time Step 

Solution 

Time Step 

Solution 

(s) 

Time (s) 

(s) 

Time (s) 

(s) 

Time (s) 

0.06 

12300 

5 

782 

Variable: 

ni 



10 

569 

15-93.4 




25 

507 

i 

- ■■■ J( 

29-172 



Solution Times on Langley CDC CYBER 173 Computer System 

* Specified relative error tolerance 0.001 
** Specified relative error tolerance 0.01 


24 















Table 5. - Effect of Modeling on Solution Times* for 
Insulated Cylinder Problem 


\ Program 

\ and 

\ Model 

SPAR (Ref. 2) 

MITAS 
(Ref 11) 

Algorithm 

Model I- 
3-D metal 
and insulation 
elements 

Model II- 
2-D metal 
elements, 3-D 
insulation 
elements 

Model III- 
2-D metal 
elements, 1-D 
insulation 
elements 

lumped 

parameter 

model 

Explicit 
(time step) 

■ 

f1 -40. ) 

77 

(4.6-40.) 

92 

(25) 

Implicit* ** 

(DT=5s) 

782 

770 

403 

387 

Implicit** 

(DT=10s) 

569 

556 

260 

1 

238 

Implicit** 

(DT=17s) 

482 

534 

216 

j.66 

Impl icit** 
(DT=25s) 

507 

Non 

Convergence 

| 

Non 

Convergence 

125 


* Time in seconds for CDC CYlJER 173 Computer 

** Crank-Nicholson and Backward Differences 
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Table 6. - Comparison of Algorithms for 'ansient Thermal 
Analysis of Titanium Multi TPS 


EXPLICIT 

crank-nickhlson 

BACKWARD 

DIFFERENCES 

I Time Step 
(S) 

Solution 
Time (s) 

Ti.'o Step 
vs) 

Solution 
Time (s) 

Time Step 
(s) 

Solution 
Time (s) 

0.03 

42,000* 

1 

1811 

1 

1772 



5 

675 

1 

5 

703 


* Estimate Based on 8400 CPU s for 400 s of Temperature History 
Solution Times for CDC CYBER 173 Computer 
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Figure 1.- Temperature history for bar example. 
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Aluminum structure 
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Figure 2.- Finite element and lumped parameter models of shuttle frame. 





Applied outer surface temperature (Node 29) 

Analysis (SPAR explicit (DT = .16s), implicit (DT = 50s)) 



o o o o o 

o o O o 

oo co ^ ea 



o o o o o o V o 

o o o o o o 

c- co m co ea 

yi ‘sjn^Bjaduiaj, 


29 


Figure 3.- Temperature history in outer structural surface of shuttle frame (Node 309). 
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Figure 4.- Finite element model of insulated cylinder. 
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Figure 5.- Heating load at outez surface of insulated cylinder. 
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Figure 6.- Effects of choice of algorithm and model changes on temperature history of 
insulated cylinder. Model I: all 3D elements. Model II: insulation- 3D, metal- 2D. 
Model Hi: insulation- ID, metal- 2D. 
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(a) Overall construction. 



(b) Representation of dimpled layer. 

Figure 7.- Titanium multiwall thermal protection system. 
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Figure 8.- Finite element model of titanium multiwall TPS. 
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Figure 10.- Finite element model of shuttle orbiter wing 
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